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Abstract. We present two numerical methods for the fully nonlinear elliptic Monge- 
Ampere equation. The first is a pseudo transient continuation method and the 
second is a pure pseudo time marching method. The methods are proved to con- 
verge for smooth solutions. We give numerical evidence that they are also able to 
capture the viscosity solution of the Monge- Ampere equation. Even in the case of 
the degenerate Monge- Ampere equation, the time marching method appears also 
to compute the viscosity solution. 



1. Introduction 

We are interested in numerical solutions of the fully nonlinear Monge- Ampere equa- 
tion 

(1.1) det D 2 u = f in f2, u = g on <9f2, 

on a convex bounded domain Q of M. n , n = 2, 3 with boundary dfl. The unknown u is 
a real valued function and /, g are given functions with / > in the non degenerate 
case and / > in the degenerate case. We will also assume that / G C(Q) and 
g G C(dQ). For / > a solution is either convex or concave. It is therefore not 



restrictive to consider only convex solutions of (1.1). 

Starting with [SJ, [IB] , interest has grown for finite element methods which are able to 
capture viscosity solutions of second order fully nonlinear equations. In the context 
of non-smooth solutions, for proven convergence results, wide stencils finite difference 
have been used for the Monge- Ampere and the Pucci equations [36] . See also [37] for a 
geometric approach. Outside of that framework, there are strong numerical evidence 
that standard finite difference methods [9] or standard finite element methods [33] 
converge as well to viscosity solutions. See also [20] and the references therein for the 
vanishing moment methodology. In this paper we give numerical evidence that C 1 



conforming approximations of a variational formulation of (1.1) converge to viscosity 
solutions and new numerical evidence that standard finite difference methods converge 
as well. This is achieved by discretizing new iterative methods we introduce. 

1.1. Contributions of this paper. We establish in the context of C l conforming 
approximations that discrete functions near a strictly convex solution are strictly 
convex. This explains why convexity did not need to be imposed explicitly in some 
previous studies. Under the assumption that a discrete strictly convex solution exi- 
sts, we prove the convergence of our iterative methods. In this paper, we prove the 



existence of a discrete strictly convex solution under the assumption that (1.1) has a 
smooth solution. In that case, Newton's method is the appropriate method for solving 

l 
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the equation. We give a proof of convergence of Newton's method when the itera- 
tes are discretized with C 1 conforming approximations. Our numerical experiments 
indicate that in the degenerate case / > 0, our iterative methods enforce convexity 
for the two dimensional problem. We prove this result at the continuous level. In 
subsequent papers, we will prove existence and uniqueness of a discrete solution in 
the degenerate case, convergence to viscosity solutions, the convexity preserving pro- 
perty of our iterative methods at the discrete level as well as introduce new iterative 
methods, with proof of convergence, which enforce convexity in three dimensions. 

The purpose of this paper is to introduce new iterative methods which are shown to 



numerically converge to viscosity solutions of (1.1) and to analyze their convergence 
for smooth solutions. Numerical results are given with the spline element method for 
which the author has extensive expertise. The finite difference results complements 
the ones obtained with the spline element method and illustrates the broad range of 
applicability of the methods presented. 

To the author's best knowledge, this is the first time the pseudo transient method 
with L = A and the truncated time marching method have been shown to numerically 
converge to viscosity solutions of the Monge- Ampere equation. 

As an initial guess for our iterative methods we take the solution of the Poisson 
equation Am = nf l l n in fi, u = g on d£l. 

1.2. Pseudo transient continuation method. The first method is a pseudo tran- 
sient continuation method. Given a fully nonlinear elliptic equation F(u) = with 
F differentiable, we consider the sequence of problems 

(vL + F'{u k ))(u k+1 - u k ) = -F(u k ), 

where L is a linear operator which can be taken as the negative of the identity or 
L = A where A is the Laplace operator and v > is a parameter. The motivation 
to consider pseudo transient continuation methods stems from the observation that 
Newton's method fails to converge to viscosity solutions. We refer to [29J for a review 
of these type of methods in the context of nonlinear equations. We presented the 
numerical results at the Thirteen International Conference on Approximation Theory 
in May 2010. In this paper, we will mainly consider the case where L = A which may 
be viewed as a preconditioner and in this case, the computations are faster compared 
to the case where L is the negative of the identity operator. Nethertheless, we will 
prove convergence of the methods at the discrete level for both cases. In the case 
of the Monge-Ampere equation, F{u) = det D 2 u — f, the solution of the problem 
is reduced to solving a sequence of elliptic equations. Recall that in that case (see 



Lemma 2.3) 



F'{u k )(u k+ i - u k ) = div {(cof D 2 u k )(Du k+1 - Du k ) = (cof D 2 u k ) : (D 2 u k+1 - D 2 u k ), 

where cof D 2 u denotes the matrix of cofactors of the Hessian and A : B denotes the 
Frobenius product of the matrices A and B. 

Given an initial guess Uq, the pseudo transient method for the Monge-Ampere equa- 
tion consists in solving the sequence of approximate problems 

(1.2) vL6 k + (cof D 2 u k ) : D 2 9 k = (/ - f k ), f k = det D 2 u k , 9 k = u k+1 - u k . 
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1.3. Pseudo time marching method. The second method is a pseudo time mar- 
ching method. Given v > 0, we consider the sequence of iterates 

(1.3) - vAu k+1 = -vAu k + F(u k ), u k+ i = g on d£l. 

This can be interpreted as an Euler discretization of the pseudo time dependent 
equation ^Jf + F(u) = 0, or as a Laplacian preconditioner of a simple pseudo time 
marching algorithm, [2B] u k+ i = u k — ^ / A~ 1 F(u k ). See also a remark in [31]. The 
simple pseudo time marching algorithm also performs well for numerical solutions in 
some cases for v sufficiently large. However the use of the Laplacian preconditioner 
besides the speed of computation, also helps select a convex solution for the two 
dimensional Monge- Ampere equation. 

The method we recommend for the degenerate Monge-Ampere equation (especially 



in three dimensions) is a truncated version of (1.3). For m = 1,2, . . ., we consider 
truncating functions Xm(%) defined by Xm(x) = —m for x < —m, Xm(%) = x for 
— m < x < m and Xm{x) = m for x > m and the sequence of problems 

(1.4) - vAuZ-i = Xm(-i/A< + F«)), t^+i = 9 on 00. 

We give numerical evidence that the limit of the discrete approximation of u™ as 



k — > oo and m — > oo approximates the viscosity solution of (1.1). Note that if F 



is smooth and the solution of (1.1) is also smooth, the right hand side of (1.4) is 



bounded and the method essentially reduces to (1.3). We have found (1.4) effective 
in the degenerate case / > but / > in Q and in three dimensions with the finite 
difference method. 

1.4. Positivity preservation of the Laplacian. We will see that the time mar- 
ching method (with preconditioner) and the pseudo transient continuation method, 
enforce Au > 0, which when combined with det D 2 u — f > gives convexity for the 
two dimensional problem. We establish this result at the continuous level. The result 
at the continuous level explains but does not prove in the degenerate case why we 
also observe numerical convexity with the spline element method. Indeed a piecewise 
convex function which is globally C 1 is convex by [30] Lemma 1. Our existence and 
uniqueness result at the discrete level cover the cases / > c > for a constant c . 
In that case we show in the paper that convexity is automatically preserved in some 



neighborhood of the solution (see Lemma 2.5). The use of the spline element method 



is also motivated by its higher order of accuracy and its robustness in some limiting 
situations [3j. 



1.5. Advantages and comparison of the two methods. The advantage of (1.4) 
over the method of Oberman [36] is that standard finite difference methods can be 
used. Potentially, the use of wide stencils can be avoided. It was pointed out in 
[22] that the simple pseudo time marching algorithm suffers from a severe time-step 
restriction, is not suitable for solutions which are not strictly convex and does not 



enforce convexity. The truncated time marching method (1.4) numerically preserves 
convexity in two dimensions and for singular solutions appears to be more accurate 
than the pseudo transient methods for small values of the mesh size. Although the 
theory of the Monge-Ampere equation has concentrated on convex solutions, one can 



equally focus on concave solutions. We found out that (1.4) is better able to capture 
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concave solutions. It is easy to implement, requiring only a Poisson solver. For 
example one can capture weak solutions of the Monge-Ampere equation by simply 



discretizing (1.4) with the standard Lagrange finite elements. The time marching 
method can also be applied to fully nonlinear equations such as the Pucci equation 
where F is not differentiate. In summary the pseudo transient methods are better 
for smooth solutions and singular solutions on a coarse mesh. Otherwise the method 
of choice is the truncated time marching method. 

The methods we propose can be used in the context of different types of discretizations 
allowing us in particular to treat more easily non-rectangular domains. The methods 
can be accelerated with fast Poisson solvers and multigrid methods. This latter 
property is even more striking for the time marching method as its implementation 
requires only having access to a multigrid Poisson solver. 



1.6. Newton's method. It is believed that some methods are able to capture the 
viscosity solutions and others cannot. This paper gives evidence that at least in two 
dimensions, by relaxation and truncation, a method which works only for solutions 
in H 2 (Q) can be expected to perform in the non-smooth case for a class of problems 
which include the Monge-Ampere equation. 

If one is only interested in smooth solutions of the Monge-Ampere equation, a num- 
ber of methods have been recently proposed. In particular, Bohmer's method JTU] 
has been recently implemented in [H] and Brenner et al [H] have used the efficient 
discontinuous Galerkin formulations. However, in all of the above cited work, New- 
ton's method is used to solve the nonlinear system of equations resulting from the 
discretization. We postulate that solving the nonlinear equations with time mar- 
ching methods captures weak solutions of the Monge-Ampere equation. In the case 
of conforming discretizations, this is evidenced by the numerical results of this paper. 
Since we are able to capture weak solutions with that technique with standard finite 
elements, the same should be true with the methods in [TTJ. Note that Newton's 
method has been used effectively in [21]. There, Newton's method is applied after 
the discretization of the equation. In our approach, the discretization process is ap- 
plied on the Newton-Kantorovich iterations. Moreover in [21] a discretization which 
is known to converge to viscosity solutions is used on parts of the domain where the 
solution is not smooth. The term uAu in the methods we propose can be seen as an 
iterative regularization term. Pryer and Lakkis [33] have shown how to use Lagrange 
elements in a nonvariational formulation for the Monge-Ampere equation but their 
method has limited applications to more general fully nonlinear equations as they 
seem able to handle only the homogeneous Pucci equation. 



1.7. Main idea. Our approach builds on the observation that globally C 1 spline 
approximations of a smooth strictly convex function remain strictly convex if the 
mesh is sufficiently fine, see Lemma 2.5 but also [3D]. This allows us to transfer well 



established results for nonlinear elliptic equations, e.g. [19] . to the context of smooth 
solutions of the elliptic Monge-Ampere equation. In this paper, we prove convergence 
results for the iterative methods introduced above from that point of view. 
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1.8. Organization of the paper. The paper is organized as follows: in the second 
section, we first prove the convergence of the methods in Holder spaces and show 
how they preserve positivity of the Laplacian and hence convexity in two dimensions. 
We then study the convergence in Sobolev spaces for smooth solutions in the third 
section. After some preliminaries, we show that a discrete strictly convex solution 



exists when (1.1) has a solution which is C 2 up to the boundary. We give error 
estimates and convergence of Newton's method. The proof of convergence of our 
iterative methods conclude the third section. The last section is devoted to numerical 
results with the spline element method and finite difference methods. It includes also 
a brief description of the spline element method. 

2. Convergence of the iterative methods 

In this section, we first describe the convergence properties of the methods in Holder 
spaces and establish that they preserve the positivity of the Laplacian. We then give 
convergence results at the discrete level for smooth solutions in Sobolev spaces. 

2.1. Convergence of pseudo transient continuation methods in Holder spaces. 



We consider a damped version of (1.2), namely 



{vL + F'(u k ))(u k+ i - u k ) = —F(u k ), 
(2-1) ^ _ 

vLQ k + (cof D 2 u k ) : D 2 6 k = -(/ - f k ), f k = det D 2 u k , 9 k = u k+1 - u k , 

T 

where r > is a damping parameter. In the numerical experiments, we used r = 1. 
We will need the following global regularity result, |40j . 

Theorem 2.1. Let Q be a uniformly convex domain in W 1 , with boundary in C 3 . 



Suppose g G C 3 (Q), inf f > 0, and f G C a for some a G (0, 1). Then (1.1) has a 
convex solution u which satisfies the a priori estimate 

\\ u \\c 2 ' a {Q.) ^ C> 

where C depends only onn,a, inf f, Q, \ \f\ \c°>fm an d IMIc 3 - 

According to [3D], all assumptions in the above theorem are sharp. We have the 
following analogue of Theorem 2.1 in [34J. 

Theorem 2.2. Let Q be a uniformly convex domain in M. n , with boundary in C 3 . 
Let < m < f < M, f G C a for some m, M > and a G (0, 1). Assume also that 
g G C 3 (Q). Then, there exists r > 1 depending on m, f , such that if u k is the sequence 
defined by ( |2.1 ), it converges in C 2 ' 13 to a solution u of (1.1) for uq sufficiently close 



to u and for every (3 < a. 

Proof. The proof follows essentially [33] • One shows by induction that there are 
constants Ci, C2 > such that 

(2.2) det D 2 u < det D 2 u k < d det D 2 u and 1 1 det D 2 u - det D 2 u k \ | c „ < C 2 , 

for r sufficiently large and v in an appropriate range. This implies that the sequence 
f k = det D 2 u k is bounded in C a and by Theorem 2.1, the sequence u k is bounded in 
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C ,a . Arzela-Ascoli's theorem is then used to prove that the sequence is precompact 



in C 2,l3 ,(3 < a. Since (1.1) has at most two solutions, PS] p. 324, by requiring u 



sufficiently close to u, we assure that the solution is locally unique. □ 

We now explain how the condition Am > is preserved in the iteration 

(2.3) uAuk+i + (cof D 2 Uk) ■ D 2 Uk+i = vAuk + (cof D 2 Uk) : D 2 Uk — det D 2 Uk + /. 

Recall that the eigenvalues are continuous functions of the entries of D 2 v as roots of 
the characteristic equation, [38] Appendix K, or [27]. Hence for 0<m<f<Mas 
in the theorem, v is strictly convex whenever \\u — i>||c a >0(n) is sufficiently small. Since 
the sequence Uk was shown to converge in C 2,l3 (Q), (3 < a, convexity is automatically 
enforced not only for u but for Uk for k sufficiently large. 

Now assume that / > and that the sequence Uk has been shown to converge to u 
in C 2,l3 (Q) for some /3 in (0, 1). From the arithmetic-geometric inequality, we have 

^!>detDV 

Again, by the continuity of the eigenvalues, Av is bounded in a neighborhood of u in 
which all Uk belong for k large enough. Choose v such that v > (n — l){Auk) n ~ l / n n 
for all k and note that the right hand of ( 2.3[ ) is equal to vAuk + (n — 1) det D 2 Uk + f ■ 



By the assumption on z/, we get vAitk+i + (cof D 2 iik) ■ D 2 Uk+i > 0. In the limit, we 
obtain uAu + (cof D 2 u) : D 2 u > 0. Since det D 2 u > by assumption, we get Au > 0. 

As for the time marching method 

— uAuk+i = —vAiik + det D 2 Uk — f, u^ + i = g on dQ, 

assume now again that / > and that the sequence Uk has been shown to converge 
to u in C 2,/3 (f2) for some (3 in (0,1). Choose v such that v > (AM fe ) n_1 /ra n . We 
have —vAuk + i = —vAu^ + det D 2 Uk — f ■ It follows from the arithmetic-geometric 
inequality that 

vAu k > — > det D 2 u k . 

n n 

and so —vAuk + det D 2 Uk < and it follows that the time marching method also 
preserves the positivity of the Laplacian. 

2.2. Convergence at the discrete level in Sobolev spaces. We first give a vari- 



ational formulation of (1.1) and its discretization by conforming finite dimensional 



spaces. Then we prove a number of lemmas which are needed for our convergence 



proofs. In particular we establish in Lemma 2J3 that discrete convex functions near 
a strictly convex solution are strictly convex. We also show the existence of a strictly 
discrete convex solution when the equation has a strictly convex smooth solution and 
give error estimates. We prove convergence of Newton's method and of our iterative 
methods. 

We use the standard notation of Sobolev spaces W k,p (Q) with norms ||.||fc iP and semi- 
norm |.|fe, p . In particular, H k (Q) = W k ' 2 (Q) and in this case, the norm and semi-norms 
will be denoted respectively by \ \-\\k and semi- norm \ .\k- We make the usual convention 
of denoting constants by C but will occasionally index some constants. For constants 
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which depend on the mesh size, we use c(h) or C(h). We make the assumption that 
the boundary of Q is polygonal and that the triangulation is shape regular in the 
sense that there is a constant C > such that any triangle K, hx/px < C, where hx 
denotes the diameter of K and px the radius of the largest ball contained in K. We 
also require the triangulation to be quasi-uniform in the sense that h/h m i n is bounded 
where h and h m i n are the maximum and minimum respectively of {hx, K G Th\. 

We choose 

(2.4) V h := S r d (T) = {se C r (n), s \ t g P d , Vt e T}. 

for a triangulation T of the domain and Vd the space of polynomials of degree less 
than or equal to d. In two dimensions, it is known that, [32], for d > 3r + 2 and 
< Z < cZ, there exists a linear quasi-interpolation operator Q mapping Li(Q) into 
the spline space S r d {T) and a constant C such that if / is in the Sobolev space 
W l+l *(n),l <p< oo 

(2.5) \\f-Qf\\k, P <Ch l+1 - k \f\ l+1 , p , 

for < k < I. If Q is convex, the constant C depends only on d, I and on the smallest 
angle 9h in T ■ In the nonconvex case, C depends only on the Lipschitz constant 
associated with the boundary of Q. 

In three dimensions, the result holds for d > 8r + l, c.f. [32]. It is also known c.f. [IT] 
that the full approximation property for spline spaces holds for certain combinations 



of d and r on special triangulations. Note that, by (2.5), 



(2.6) \\Qu\\ 2 , p < C\\u\\ 2:P , ueW^ip), 

for all p. We assume that r > 1 and that the following inverse inequality which holds 
for finite element spaces, c.f. Theorem 4.5.11 of [12], also holds for the spline spaces 



(2.7) \\u\\ 2 , p < Ch l - 2+ ™ n ^-^\\u\\ l>q ,\/u G V\ 

for < Z < 2, 1 < p, q < oo. The local estimates may be viewed as a consequence of 
the assumption of uniform triangulation and of Markov inequality, [32] p. 2. Passing 
from local estimates to global estimates can be done as in |31j . 



The variational formulation of (1.1) is given by: find u G W 2 ' n (Q), u = g on dQ such 
that 

(2.8) - - I (cof D 2 u)Du -Dw dx = [ fw dx, Vu> G W 2 ' n {Vt) n H 1 ^). 
n Jn Jn 

To see that the left hand side is bounded for u G W 2,n (Q), notice that for n = 2 



(cof D 2 u)Du ■ Dw dx 



n 



< C\\D 2 u\\ ,2\\Du\\ 0A \\Dw\\ 0j 4. 



Next for u G H 2 (Vt), du/dxi G i/ 1 (fi), i = 1, . . . , n and by the embedding of if x (f2) in 
L g (Q) for q > 1 when n = 2, the right hand side above is bounded by C\ \D 2 u\ \L 2 (n)\ \ u \ \h 2 (q) 
\\ w \\H 2 (n)- in three dimensions, the term cof D 2 u involves the product of two second 
order derivatives. We have by Holder's inequality and the embedding of if 1 (Q) in 
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L q (Q) for 1 < q < 6 when n = 3 

f d 2 u d 2 u du du 
dx 2 dz 2 dx dy 



dx 



< 



d 2 w 
dx 2 



1 0,3 | 



d 2 u 
dx 2 



0,3 



d 2 u 
dx 2 



0,6 



d 2 w 
dx 2 



0,6 



We conclude that for n = 3, 

/ (cof D 2 u)Du ■ Dw dx 
Jn 

We note that for n = 2, 3, we may write 



(2.9) 



(coi D 2 u)Du ■ Dw dx 



< C\\u\\l 3 \\u\\ 2 \\w\\2. 



< ciMI^IMhlMk 



Put V = W 2 ' n {Q) and V = W 2 > n {Q) n H^{Q). Let V h be a conforming finite di- 
mensional subspace of W 2,n (Q), Vq be a conforming finite dimensional subspace of 
W 2,n (Q) PI /fg(^)- Furthermore let ^ be the interpolant in ^ of a continuous ex- 
tension of g. We have the following conforming discretization of (2.8): find G V h , 
u>h — 9h on dVL such that 



(2.10) 



— / (cof D 2 Uh)Dui l ■ Dwh dx = [ fwh dx, Wwh & Vq. 
n Jn Jn 



(2.11) 

And for F(u) 



We now prove a number of preliminary results. The following lemma is elementary. 
Lemma 2.3. We have 

det D 2 u = -(cofD 2 u) : D 2 u = - div ((cofD 2 u)Du). 
n n 

det D 2 u we have 
F'(u)(w) = (cofD 2 u) : D 2 w = div ((cofD 2 u)Dw), 
for u, w sufficiently smooth. 

Proof. Note that for any n x n matrix A, det A = l/n(cof A) : A, where cof A is the 
matrix of cofactors of A and for two nxn matrices M, N, M : N = Y^j=i M^iVy. For 
any sufficiently smooth matrix field A and vector field v, div A T v = (div A)-v+A : Dv. 
Here the divergence of a matrix field is the divergence operator applied row-wise. If 
we put v = Du, then det_D 2 w = l/n(cof D 2 u) : D 2 u = 1/n (cof Dv) : Dv and 
div(cof Dv) T v = div(cof-Du) • v + (cofDv) : Dv. But div cof Dv = 0, c.f. for example 
[TS] p. 440. Hence since D 2 u and cof D 2 u are symmetric matrices (2.11) follows. The 
assertion about the Frechet derivative of F also follows from these identities. □ 

Recall that W 2,P (Q) is continuously embedded in H 2 (Q) for p > 2 and hence for 
ueW 2 ' n (n),n>2 



(2.12) 



\u\ 



< C\\u\\ 2 ,n- 



Lemma 2.4. Letv,w G W 2 ' n (n),n = 2,3 and $ G n H 2 {9), th 



en 



(det D 2 v - det D 2 w)t/> 



((cof[(l -t)D 2 w + tD 2 v)) 
(Dv-Dw)) -DipdxXdt, 



and 
(2.13) 



/ (det D 2 v - det D 2 w)tp dx 
Jn 



Moreover, i/v,ty G W 2 ' n (Q) n W 2 '°°(fi) 



(2.14) 



(det -D 2 -y - det D 2 w)-0 dx 



< Co(||v||2,oo + lk|| 2 ,oo) n \\v- 



Proof. We first recall the Mean Value Theorem. Let E and F be Banach spaces and 
and let us denote by L(E,F) the space of continuous linear mappings from E to F. 
Let also X be an open subset of E and let F : X — > F be a differentiate map. If 
F' : X — > L(E, F) is continuous, F is said to be of class C 1 and for all a, x G X, we 
have 

F(x)=F(a)+ / F'[(l-t)a + tx](x-a)rft. 
Jo 

Next, let F : C°°(fi) C°°(fi) denote the mapping v H> det-DV Then F is 
different iable with 

= (cofD 2 w) : D 2 w = div ((cof D 2 u)Dv) . 
Since v H- -F'[t>] is linear, F is of class C 1 and by the Mean Value Theorem 

F(v)-F(w)= I div ((cof (1 -t)D 2 w + tD 2 v){Dv - Dw)) dt. 
Jo 

Next, let ip £ T>(Q). By Fubini's theorem, 

I (det D 2 v — det D 2 w)ip dx = I if div ((cof (1 - t)D 2 w + tD 2 v)(Dv - Dw))ip dx\ dt 
Jn Jo lid J 

= ~L{I (^ o{ ^~^ D2w + tD2v ^ Dv ~ Dw ^)' D ^ dx ^ dt - 

As in the proof of ( |2.9 ), we conclude that (2.13) holds for u,ip smooth, ip = on 
dQ. Then by the density of V(Q) in #<J(fi) it also holds for v,w G C°°(fi) and 
V> G i^o(^) ^ H 2 (Q). Finally by the density of T>(Q) in the Sobolev spaces, it is also 
valid for v,w G W /2 ' n (f2). Inequality (2.14) is proved similarly. □ 

The following lemma is fundamental. 

Lemma 2.5. Let u be a C 2 (Q) strictly convex function. Then there exists 5 > such 
that for ||f — w||2,oo <S, v is strictly convex. Moreover if u G W 2,OQ (fl) nH l+1 (Q), 1 + 
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n/2 < I < d, there exists 5(h) > such that for Vh G V h , \ \vh — u\\i < 5(h) and h 
sufficiently small, Vh is strictly convex and 



(2.15) 



\(coiD 2 v h (x))Dw(x)} -Dw(x) dx < M\\w\\l,w G H^(Q), 



for constants m, M > which are independent of h. It follows that for h sufficiently 
small, r > 1, Qu is strictly convex. 

Proof. Let X\(D 2 u(x)) and \2(D 2 u(x)) be the smallest and largest eigenvalues of 
D 2 u(x) respectively. Since u is strictly convex, we have 

2m < AxpVa;)) < X 2 (D 2 u(x)) < M/2, for m,M > 0. 

Since the eigenvalues are continuous functions of the entries of the Hessian, as roots of 
the characteristic equation, for all e > 0, there exists 5 > such that \ \v — w||c 2 (n) — $ 
implies \Xi(D 2 v(x)) - \i(D 2 u(x))\ < e,i = l,2,Vx G Q. Since \i(D 2 u(x)) > 2m, Vx G 
ft, for e = m, we get \\(D 2 v(x)) > m,Vx G Q. We conclude that for \ \v — M||c 2 (n) < ^> 
Xi(D 2 v (x)) > m,Vx G Q. Similarly, \2(D 2 v(x)) < M,Vx G Q. This implies that v is 
strictly convex for ||i> — w||c 2 (n) < 

Since \\(D 2 v(x)) and \2(D 2 v (x)) are the minimum and maximum of the Rayleigh 
quotient (cof D 2 v(x)z) ■ z/\\z\\ 2 , where \\z\\ denotes the standard Euclidean norm in 



we have 



m\\z\\ 2 < [(coiD 2 v(x))z\ ■ z < M\\z\\ 2 ,z G 



This implies 

(2.16) m\\w\\l < / \(coiD 2 v(x))Dw(x)} ■ Dw(x) dx < M\\w\\\,w G H^(Q), 

Jn 

where we used the equivalence of the seminorm | |i and norm || ||i on Hq(Q) Next, on 
each element K, there is 5(K) > such that for \\vh — u\\c^(k) < (2.16) holds 

with the domain Q replaced by K. Let then 5(h) be the minimum of {5(K),K in 
%}■ Since by dZ5b and for I > 2 



\Vh - u\\ 2 ,oo <\\Vh — Qu\\ 2 

C 



< 
< 
< 



h 1+ % 

c 



V \\Qu - u\ 
v h - Qu\\i + Oh 1 ' 1 



2,oo 



h 



C 



\Vh - M l + 



c 



h 



1+2 



U\l+l,oo 

u - Qu\\i + Ch l ~ 1 \u\i + i 



nWvh-uWt + Ch 1 1 r i\u\i+i + Ch l l \u\ i+1) 



i-ii 



h 1+ 



\C 2 (K) 



< 



we see that if \\Vh — «||i < h 1+n / 2 5(h) / (3C) and h sufficiently small, \\vh — u 
\\ v h — u\\2,oo < 5(h) and Vh is piecewise convex. Since r > 1, Vh is of class C 1 and 
piecewise convex that is, Vh is convex by [30], Lemma 1. 

Since by (2.9), J n [(cof D 2 Vh(x))Dw(x)] ■ Dw(x) dx is finite for w G H^(fl) and Vh G 
V h c W 2,n( Q ^ we conc i u d e that plB] holds. 

□ 
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2.2.1. Existence, uniqueness and error estimates for C 1 conforming approximations. 



We show in this section that if (1.1) has a smooth strictly convex solution, then the 
discrete equations (2.10) also have a unique strictly convex discrete solution. The 



existence and uniqueness of the solution of (2.10) was also addressed by Bohmer 



[10] from a different perspective. Bohmer's method requires stable splitting of spline 
spaces [15] and the required modification of the Argyris space has been implemented 
only recently [14] . We adapt and simplify the methods of [21] which were given for 
the vanishing moment formulation. In addition to existence and uniqueness, this also 
gives error estimates. See also [33] for another discussion of conforming approximation 
for homogeneous boundary conditions. The latter study has been completed after 
the results of this work were announced. We note that one may be able to prove the 
existence of a solution at the discrete level following [39] as indicated in the proof of 
Theorem 5.8 of |19|. 



Lemma 2.6. Assume u G W +1, °°(tt), 3 < I < d is a strictly convex function, that 
Q is convex with a Lipschitz continuous boundary and that the spaces V h = S r d (T) 
have the optimal approximation property (2.5) and satisfy the inverse estimates (2.7). 



Then the problem (2.10) has a unique strictly convex solution Uh and we have the error 
estimates 

\\u - u h \\ 2 < Ch l ~ l \u\i + x, \\u - Uh\\i < Ch l \u\i + i 
\\ u — Uh\\L 2 = 0(h l+1 ) for h sufficiently small. 

Proof. We consider the linear problem: find v G H^(Q) such that 

(2.17) / (cof D 2 u)Dv ■ Dw dx = <pw dx,\/w G H^{Q), 

Jn Jn 

for G L 2 (Q). Since u G W l+1 '°°(Q),3 < I < d is strictly convex, 3 m, M > such 
that 

m\\w\\\< i \(coi D 2 u(x))Dw(x)} ■ Dw(x) dx < M\\w\\\,w G H^(Q). 
Jn 



The existence and uniqueness of a solution to (2.17) follows immediately from Lax- 



Milgram lemma. Similarly, there exists a unique solution to the problem: Find 
v h G Vq such that 

(2.18) / (cof D 2 u)Dv h ■ Dw h dx = <j>w h dx,Vw h G V h . 
Jn Jn 

We note that the constant C above depends on u and that since Q is assumed convex, 
v G H 2 {yi) by elliptic regularity. 

We define a bilinear form on Hq(Q) x Hq(Q) by 

(2.19) B[v,w] = / {cof D 2 u)Dv ■ Dw dx, 

Jn 

and for a given Vh G V h , Vh = gn on <9f2, define T{vh) as the unique solution of 

(2.20) B[v h - T{v h )M = - [ det D 2 v h ip h dx + f f^ h dx, G V h . 

Jn Jn 
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Since Vh — T(vh) G Vq, T(vh) G V h , T(vh) = Qh on dfl. A fixed point of the nonlinear 
operator T corresponds to a solution of (2.10) and conversely if Vh is a solution of 
(2.10), then Vh is a fixed point of T. We will show that T has a unique fixed point in 
a neighborhood of Q(u). Put 

£ ft (p) = {v h G V h ,v h = g h ondtt, \\v h - Qu\\i < p}. 

We first show that 

(2.21) \\Qu - T(Qu)\\i < d^ll^ll^H^, 

then we show there exists < po depending on h such that T is a contraction mapping 
in the ball Bh(p ) with a contraction factor 1/2. We conclude by applying the Brouwer 
fixed point theorem in a suitable ball. 

Put Wh = Qu — T(Qu). Then Wh G Hq(Q) and using det D 2 u = f, 

B[wh,Wh]= / ( det D 2 u — det D 2 Qu)wh dx. 
Jn 



Then by Lemma 2.4, the coercivity of B on Hq(Q), ( 2.6[ ) and (2.5), 



Wh\\i < C \\ U \\2oo\\ U - Q U \\l\\ W h\\l < 



i||„.ltn-l 
2,oo 



from which ( 2.21[ ) follows. 

For Vh,Wh G B h (p ), with p yet to be determined, and ^ G V Q h , 
B[T(v h ) - T(w h ),'ip h ) = B[T(v h ) - v h , ip h ) + B[v h - w h , ip h ) + B[w h - T(w h ), ip h ) 

(det D 2 v h — det D 2 Wh)il>h dx 



+ / [(cof) D 2 u{Dv h - Dw h ) ■ Di) h dx. 
Jn 



Using 



/ cof D 2 u(Dv h — Dwh)Dijj h dx = / / [(cof D 2 u)(Dv h — Dwh)] ■ Dijj h dx dt, 
Jn Jo Jn 

Lemma 2.4, and with the observation that by the inverse inequality, Vh,Wh G H 2 (Q), 
we have 

B[T(v h ) - T(w h ),ij h ] = J | J[{^D 2 u ~ cof ((! " + ^V)) 

(Dfft - Dn)] • D^fr <ix| dt 
= L {I ^ CoW2u ~~ coiL>2 Q u )( Dv h - Dw h )] ■ Dij h dx^ dt 
-J {J {coi[{{l-t){D 2 w h -D 2 Qu)+t{D 2 v h -D 2 Qu)) 
(Dvh — Dwh)} ■ Dip h dx > dt. 
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We conclude using ip h = T{yh) — T(wh), the coercivity of B on Hq(Q), and (2.9), that 



\\T(v h )-T(w h )\\i < CWu-QuWr^Wvh-WhMiPh 



+ C(\\w h - Qu\\™- L + \ \v h - Q 



u 



in— IN 
\2,n J 



\Vh - W h \\ 2 \\lph\\2- 



For n = 2, by the inverse inequality (2.7) and (2.5), 



Po 



\\T(v h ) - T(w h )\\i < (C 2 ^->| m)00 + C 3 ^)\\v h - lOftlU- 

We require C , 2^~ 1 Mh-i,oo < 1/4 and we choose p Q such that, C 3 p /h 3 < 1/4 for 
example po — ^ 3 / (4C3). It follows that T is a contraction mapping in the ball Bh(p ) 
with a contraction factor 1/2. We now consider the case n = 3. We have by (2.7) 

I|tK)-tWII 1 <c^ ( '- 1) Ki. 1 



- Wh\\i 

|2 



+ C(||u; ft - Qu\\ 23 + \\v h - Qu\\ 23 ) 



\Vh - w h \\i 
h 2 



Again, by (2.7), 



We conclude that for n = 3, 



C 

\ V \\2,3 < 7T 



\v\\o < 



C 



\v\\i. 



\T(v h )-T(w h )\\ 1 <{C 4 h« l - 1) 



P 2 

^oc + ^pjlk-^lli- 



Now, we require C A h 2<yl l ^\u\ 2 +loo < 1/4 and choose p such that C 5 pl/h 5 < 1/4. 
Finally, note that with p sufficiently small, in B^po), 

\\Tvh-Qu\lx < WQu-TiQu^ + WTQu-TvhWt < C 1 h l \\u\\^\u\ l+1 + \\ Vh ~^ u W\ 

Put px = 2Cih l \\u\\%^\u\i+i. Since I > 3 and h sufficiently small, p\ < po, and T 
maps Bh(pi) into itself. 

We conclude by the Brouwer fixed point theorem that T has a unique fixed point Uh 
in Bh(pi) to which the iterates = T(u^) converge. Moreover 



\u - u h \\i < \\u - Qu\\x + \\Qu - U h \ 



\u - Qu\\i + \\Qu - T{u h )\\-i 



< Ch l \u\ l+1 + Pl < {C + 2C x \\u\\r^)ti\u\i+i 



and 



\u - u h \\ 2 < \\u — Qu\\i + \\Qu - Uh\\2 = \\u - Qu\\ 2 + \\Qu - Tuh\\i 



Pi 



<Ch l - L \u\ l+1 + ^<Cti- l \u\ l+1 , 
n 

using again an inverse estimate. 

Using a duality argument, we prove the L 2 error estimate. Recall that the domain is 
convex and let w G H 2 (Q) be the solution of the problem 

div(cof D 2 u)Dw = u — Uh, in f2, w = on 

By elliptic regularity, 



(2.22) 



Mb < C\\u — Uh\ 
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\u - u h 



(2.23) 



u — u h ) div(cof D 2 u)Dw dx 

I [(cof D 2 u)Dw] ■ D(u - u h ) dx 
Jn 

/ [(cof D 2 u)(Dw - DQw)} ■ D{u - u h ) dx 
Jn 

+ / [(cof D 2 u) DQw] ■ D[u - u h ) dx. 
Jn 



Recall that for Vh G Vq 



' det D 2 uv h dx = / [(cof D 2 u)Du] ■ Dvh dx = I fvhdx, 
n Jn Jn 



det D 2 UhVh dx 



Hence by Lemma 2.4 

= / det D 2 uVh dx — I det D 2 UhVh dx 
Jn Jn 



■ Dvh dx 


= / fv h dx. 




Jn 


[(cof(l- 


- t)D 2 u h + tD 2 u) 


(Du- 


Du h )} ■ Dv h <is| dt, 


[(cof(l- 


- t)(D 2 u h - D 2 u) + D 2 u) 


(Du- 


Duh)} ■ Dvh c?x| dt. 



Subtracting from (2.23), using Vh = Qw, we obtain 

||w - Uh\\l < C\\w - Qw\\i\\u - Uh\\i + +C\\u - w^ll^n 1 !!"" _ "ftlbllQwIh- 

Since \\w — Qw\\\ < Ch\\w\\2, |Mk — ll^lb) HQ^Ih < 1Mb, by (2.22), 

(2.24) | \u - u h \ | < Ch l+l \u\ l+l +C\\u- UhW^h 1 - 1 \u\ l+1 , 

which gives when n = 2 since I > 3, ||u — w^Ho — 0(h l+1 ). 

For n = 3, as in the proof of the H 1 error estimate using p 1; we have 



C 



\U - U h \\ 2> 3 <\\U~ Qu\\2,3 +\\Qu- U h \\ 2) 3 < h \u\l +1)3 + -j\\Qu - Tu h \\l 



l+l- 



Again, here for h sufficiently small, we get from (2.24) that the term h l+1 dominates 
for n = 3. □ 

Assume that Uh, the unique strictly convex solution of the discrete Monge-Ampere 
equation, satisfies Uh\\i < <W2, where 5h > has the property that if \\vh— u\\i < 
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Sh,Vh G V , then Vh is strictly convex. This is possible by Lemmas 2.5 and 2.6 under 
the assumption u G W l+1 '°°(Q), 1 + n/2 < I < d. We put 

X h = {v h G V h ,v h = g h ondQ,\\v h - < 5 h /2}, 



so that for Vh G X h , Vh is strictly convex. Moreover ( |2.15 ) holds. 



We consider the mapping Fh : X h — > ( Vq 1 )', where (Vq 1 )' is the dual space of V h and V h 
is equipped with ||.|| 2 , defined by (F h (v h ) , i[> h ) = f n det D 2 v h ip h dx, v h G V h ,if) h G 
and recall that / > is continuous. Since ft is bounded, by L? duality X h C (Vq)'. 
We use the notation ||.|| for the operator norm of an element of a dual space. With 



these notation, (2.8) can be written det D 2 u = f in V ' and fl2.10p can be written 
F h (u h ) = / in (Vjf. We have 

Lemma 2.7. Discrete coercivity: \\Ffi(vh)(p)\\ > c/i||p||i ; for all p G V h and Vh G X h 
for a constant c > 0. 

Generalized Lipschitz continuity: \\F' h (vh) — F' h (wh)\\ < c(h)\\vh — Wh\\ 2 n , Vh, Wh G X h 
where c(h) is independent of h for n = 2 and c(h) = Ch~ 3 / 2 (5h + \ \uh\\i) for n = 3. 
The constant c(h) depends on 5h and u^. 

Proof. Note that 

(2.25) m\\p\\\ < (F'(v h )(p),p) < M\\p\\lp G ^(O), 

v h G X h and ||J^(t;)(p)|| = sup^ | (F h (v)(p), ^) |/| H | 2 > |<F;»(p),p>|/|H| 2 > 



m(|b||i/|b|| 2 )|b||i. By (23, IWI1/IWI2 > ch for all p G V h . This proves that 

— c ^lblli)P Tq for a constant c > 0. 
For u A , G X h , *0 G X' 1 , 77 G V^, we have 

(F , h (v h )(ip),rj)-(F' h (w h )(1>) 1 ri)= [ (div(cof D 2 v h )DiP)r] dx - [ (div(cof D 2 w h )DiP)r)dx 

Jn Jn 



[(cof D 2 v h )Dip} ■Dr l dx + / [(cof D 2 w h )Dip} ■ Dr]dx 

Jn 

I [(cof D 2 w h - cof lPv h )Di/)] ■ Dr]dx 
Jn 



For n = 2 cof Dwh — cof D vh = cof D (wh — Vh) and as for (2.9), we have 



IKK)(^)-F;K)^)||<c|K-^|| 2 ||vii 2 
H^K)-^K)||<qK-^|| 2 . 

For n = 3, the integral J n [(cof D 2 Wh — cof D 2 v h) Dip]- Dr/dx is the sum of terms similar 
to 



d 2 w h d 2 w h ( d 2 w h 



dy 2 dz 2 



dydz 



d 2 w h d 2 w h ( d 2 w h 



dy 2 dz 2 



dydz 



dip dr] 
dx dx 



dx. 
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By the mean value theorem and as with (|2.9|), and using inverse estimates, we obtain 

I2IMI2 



\(K( v h)(i>),v) - (Ki w h)W,v)\ < C(\\w h \\ 2i3 + \\v h \\ 2!3 )\\v h -w h \\ 2t3 

\\K(vh)W - ^HXVOII < c(||tw fc || 2 , 3 + ikMlk - ™h\\ 2 , 3 

\\ F h( V h) - K( W h)\\ < C(\\w h \\ 2t3 + \\v h \\ 2i3 )\\v h - W h \\ 2t3 

3 

< Ch~?(5 h + - w h \\ 2>3 , 

which proves the result. 



□ 



2.2.2. Convergence of Newton's method for C 1 conforming approximations. When 
(1.1) has a smooth solution, the appropriate way to solve the nonlinear equations 
is to use Newton's method. We establish here that the discretization of the iterates 
converge to the unique local solution of (2.10). We have 



Theorem 2.8. Under the assumptions of Theorem 2.6, the sequence defined by 

F' h (u k )(u k+1 - u k ) = -(F h (u k ) - /), 



with a suitable initial guess converges to the unique convex solution Uh of (2.10) with 
a quadratic convergence rate. 



Proof. Since > c/i||p||i,p G Vq and u k+ i — u k e Vq, the Newton's iterates 

are well defined. We first show that for u k G X h , 

(2.26) \\F h (u k+1 ) -f\\< c(h)\\F h (u k ) - /|| 2 . 

Then we establish that if u k G X h u k+ \ G X h 5(h) sufficiently small. Finally we prove 
the convergence rate. 

Put p k = u k+ \ — u k and let iph G Vq. We have by the definition of Newton's method, 
Fh(uh), by Fubini's theorem and (2.9), 

(F h (u k+1 ) - f,if) h ) = (F h (u k+ i) - F h (u k ) + F h (u k ) - f,ip h ) 

= ( F' h (u k + tp k )p k dt + F h (u k ) - /, ipt) 
Jo 

= (F h (u k ) - f + F' h (u k )p k + / F' h (u k + tp k )p k - F' h (u k )p k dt, ip h ) 

Jo 



< 



\(F' h (u k + tp k )p k - F' h (u k )p k dt,i)) < c(h)\\p k \\ 2 \\p k \\ 2 \\il) h \ 



which implies \\F h (u k+ i) — f\\ < c(h)\\p k \\ 2n \\p k \\ 2 . Since u k G X h is strictly con- 



/II 



\Fi(u k )( Pk )\\ > 



vex and F' h {u k )(p k ) = -(F h (u k ) - f), we have ||i^(«fc) 

c/i||p fc ||i- And by ( J2?f| ), |k|| 2 , n < ch l ' n/2 \ \v h \ \ 2 and \\v h \\ 2 < c /h\\vh \\i. So ||^||2,n < 
ch' n / 2 \\v h \\ x . That is, |bfe|| 2jn ||pfe|| 2 < Ch-^^WpkWl. Hence ( ggg holds. Next, we 
prove that if u k G X h , then u k+ \ G X . We have 

F' h (u k )(u k+1 - u h ) + F' h (u k )(u h - Uk) = -{F h (u k ) - f). 

Hence, using F h (uh) = f and the Mean Value Theorem, 

K( u k)( u k+i ~ u h ) = -F' h {u k )(u h - u k ) - F' h [u h + 6(u k - u h ))(u k - u h ) 

= [F h {u k )-F h {u h + e(uk-u h ))}(u k -Uh),9e [0,1]. 
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We conclude by Lemma 2.7 



ch\\u k+1 - Uh\\i < c(h)\\u k -u h - 6(u k - u h )\\ 2 J\u k - u h \\ 2 < c(h)h n/2 1 \\u k - u h \\\. 
Hence 

IK+i - UhWx < c{h)\\u k - u h \\\ < c{h)5 h \\u k - u h \\ x . 

This gives the convergence rate and shows that, with 5(h) sufficiently small, u k+ i G 
X h when u k G X h . 

We now show that u k converges to Uh- Let us assume that uq is chosen so that 



no G X h and u k G X h for all k. From (2.26), we obtain 
\\F h {u k ) - /| | < c(h) ■ c(h) 2 ■ c(hf ■ ■ ■ cihf- 1 \\F h (u ) -f\f 



c(h)^\\F h (u ) - f\f = c(h)- 1 c(h)\\F h (u ) - /|| = c^rV 



where q = c(h)\\F h (u ) - /||. 

We now further make the assumption that uq is chosen so that q < 1. Then, by the 
Mean Value Theorem, \\F h (u k ) - f\ \ = \ \F h (u k ) - F h (u h )\\ = F' h (v h )(u k -u h ) for some 



Vh in X th . We then have by Lemma 2.7 



X 1 y \ 1 0& 

\u k - u h \\ x < —\\F h (u k ) - f\\ < —c(hy q 



which proves the convergence. 



□ 



2.2.3. Convergence of the pseudo transient methods for C 1 conforming approxima- 
tions. We are now in position to prove the convergence of the iterative methods 



(1.2). The discretization of (1.2) depends on the choice of L: Given v > and a suit- 



able initial guess, find u k+ i G V h ,u k+ i = goiadQ such that we have for all iph G Vq 1 , 
when L is the Laplace operator 

(2.27) 



-v / (Du k+X - Du k ) ■ Dtp h dx + (F' h (u k )(u k+1 - u k ),tp h ) 
Jn 



(-(F h (u k ) - f),ip h ), 



and when L is the negative of the identity, 



(2.28) 



-v I (Wfe+i - u k )ip h dx + (F' h (u k )(u k+1 - u k ),ip h ) = (-(F h (u k ) - f),ip h )- 
Jn 

Above, we have made the abuse of notation of denoting by both u k the solution of 
the iterative methods at both the continuous and discrete level. In the remainder of 
this paper, only discrete solutions are considered. This alleviates the notation. 

Theorem 2.9. Let Q be convex with a Lipschitz continuous boundary and assume that 
the spaces V = S r d (l~) have the optimal approximation property (2.5) and satisfy the 



inverse estimates ( |2.7 ). A sequence defined by either (2.27) or (2J28J) with a suitable 
initial guess and a suitable value of v converges to the unique strictly convex solution 



of (2.10) for h sufficiently small. Moreover the convergence rate is linear. 
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Proof. Define Mi : V h ->■ (V^)', i = 1, 2 for u, ^ G by 

(Mi(v),ip h ) = / Dv-D^ h dx, (M 2 (v),ip h }= / uVfcda;. 

We note that 

(2.29) \\Mi{v)\\ < \\v\\ 1 ,veV h ,i = l,2. 

Next, for p G Vq 1 , we have by the inverse inequality 

-^i(p)W + W(p)WI 



- su P^o 
I — z> 

> 



HV'fclla 

vM x (p)(p) + F h (v h )(p)(p)\ 



INI2 

+ , IpI 



> z/|p|i-j— — : — h m 
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Similarly 



WW 2 (p) + F((^)(p)||>H|p||oS 



We therefore have 

(2-30) Mix < ^^J l - vMi(p) + F' h (v h )(p)\\,p G ttf.i = 1,2. 

We can now determine under which conditions when u k G X ft we have Mfc+i G X ft as 
well. Using Fh(uh) = f and the Mean Value Theorem, 

-vMi(u k+1 - u h ) + F' h (u k ){u k+ i - u h ) = -vMi{u k - u h ) + F' h (u k ){u k - u h ) 

- F' h (u h + 9(u k - u h ))(u k - u h ) 
= \F' h (u k ) - F' h (u h + 9(u k - u h ))](u k - u h ) 
- uMi(u k -u h ),e G [0,1]. 



Using (2.30) and (2.29), the generalized Lipschitz continuity property of F' h and in- 
verse inequalities, we get 

(uQh 1 + C 3 h)\\u k+1 - u h \\i < v\\u k+1 - u h \\i + c(h)\\u k+1 - Uh|| 2l „||ufc+i - u h \\ 2 

< v\\u k -u h \\i +c(ti)\\u k -u k \\l 
<(u + c(h)S(h))\\u k - u h \\i. 

We therefore have 

. M u + c(h)5(h) .. 

(2.31) i< Mfc-^fr i- 

z/Cj/i* + 63/1 

We note that (1/ + c(h)5(h))/(uC i h i + C 3 h) < 1 is equivalent to v(l - QH 1 ) < C 3 h - 
c(h)5(h). This shows that if we choose h sufficiently small so that 1 — Cih 1 > 0, and 
then 5(h) sufficiently small so that C 3 h — c(h)5(h) > and then v sufficiently small 
and possibly 5(h) smaller, Uk+i G X h when u k G X h . 
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We now assume that Uq is chosen in X h . We have for % = 1, 2 

(F h (u k +i) - f,4>h) = (F h (u k +i) - F h (u k ) + F h (u k ) - f,ip h ) 

= (F h {u k+1 ) - F h (u k ), ip h ) - (vMi{u k+ i - u k ), if) h ) 

- {F' h (u k ){u k+ i - u k ), ip h ) 



= ( / [F' h (u k + t(u k+ i - u k )) - F' h (u k )](u k+ i - u k ) dt 
Jo 

- (uMi(u k+1 - u k ),il) h ). 

We conclude that 

\\F h (u k+ i) - f\\ < c(h)\\u k+1 - u k \\ 2jn \\uk+i - u k \\ 2 + v\\Mi(u k+ i - u k )\\. 

By ( |2~7l ), \\u k+1 - Mfc|| 2 ,nlK+i - u k \\ 2 < c(h)\\u k+1 - u k \\l. Moreover by d2.29[ ), 
\\-Mi(u k+1 — u k )\ \ < \ \u k+ i — u k \\i, i = 1,2. We therefore have 

\\F h (u k+1 ) - f\\ < c{h)\\u k+1 - u k \\\ + v\\u k+ i - u k \\ x . 

Finally, by ( 2.30[ ) and the definition of the iterative methods (2.27) and (2.28), \\u k — 



u h \\\ < u Cihi+c 3 h W Fh ( Uk ) ~ = !> 2 - We conclude that 

c(h) 



\F h (u k+1 ) - f\\ < 



\F h (u k ) - f\\ 2 + 



vdW + C 3 h 



\F h (u k )-f\\ 



Put a(is } h,n) = c(h) / (vdh 1 + C 3 h) 2 and 0(v,h,n) = v / '(udh* + C 3 h) . We note that 
(3(v, h, n) = v I [yCih 1 + C 3 h) < 1 is equivalent to v{2 — Cih % ) < C 3 h. We assume that 
h is chosen sufficiently small so that (2 — Cih 1 ) > and v is chosen sufficiently small 
so that (3(u, h,n) < 1/2. Let q = \\Fh(uo) — f\\ and assume that uo is chosen so that 
a(u, h,n)q < 1/2 and q < 1. We then have 

s = a{y, h, n)q + f3{v, h, n) < 1. 

It follows that 

\\F h (ui) -f\\< a(u, h,n)\\F h (u ) - f\\ 2 + (3(v,h,n)\\F h (u ) -f\\ = sq, 
and since s < 1, 

\\F h (u 2 )-f\\ < \\F h ( Ul ) - f\\{a{u,h,n)\\F h ( Ul ) - f\\ + P(v,h,n)) 
< \\F h ( Ul ) - f\\s < s 2 q. 

We conclude that \\F h {u k ) — f\ \ < s k q. Using F h {uh) = f, the Mean Value Theorem 
and the discrete coercivity property of F' h , we have 

1 



\u k - u h \\i < 



ch 



\F h (u k )-f\\<^s k , 



from which the convergence follows. The convergence rate is given by (2.31). 



□ 



Remark 2.10. The analysis above does not indicate whether (2.27) should be pre- 
ferred over (2.28). We view (2.27) as a preconditioned version of (2.28). Moreover, 
the numerical results indicate that the use of the Laplacian preconditioner improves 
the convexity property of the numerical solution. 
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2.2.4. Convergence of the time marching methods for C 1 conforming approximations. 
We now turn to the proof of one of the main results of this paper, the convergence 
analysis of the iterative method (1.3) for the Monge- Ampere equation. 

By Poincare's inequality, \w\\ < \\w\\i < C^\w\i,w G Hq(Q). Let •jh = bh/C^. We 
now assume that Uh the unique strictly convex solution of the discrete Monge- Ampere 
equation satisfies \u — Uh\\ < 7/i/2 and we now define X h as 

X h = {v h G V h ,v h = g h ondVt,\v h -u h \i <7^/2}. 



Then for v h G X h , by Lemma 2.5, m|w|f < J n [(cof D 2 v h )Dw] ■ Dw dx < C 2 M\w\\, 
w G Hq(P). Without loss of generality, relabeling C\M as M, we will assume that 
for Vh G X h , 

(2.32) m\w\\ < / [{cof D 2 v h )Dw] ■ Dw dx < M\w\\,w G H^Q). 

Jn 

Let v = (M + m)/2 and define a mapping T\ : X h — > (V h )' by 

/ Dv h -Dip h dx+- I (detD 2 v h -f)i) h dx,v h eX h ,'ip h eV c 
Jn v Jn 

The following lemma will make it possible to show that T\ is a strict contraction. 
Lemma 2.11. For v h G X h , 

Ml _ \\Tj(v h )(M\ , \T[(vMM h )\ ^ M-m 

N J iWN* - su P^ h ev h ,^o i^-J^ - su Pi> h ev h ,4> h ^o |^|2 - M + m' 



(2.33) {T x (v h )^ h ) 



h 

• 



Proo/. Let a = rop^^.^ ^^f^ 1 . We have 

(2.34) |T{K)WWI < e Vf. 

Since for /i h G V^, = sup^ e ^ ^ 1 ) ) C 7 ?^ ) I / 1 ^ 1 1 ? we obtain 



\ T l{ v h)\\* - su Pn h ,v h £V£,n h ,r, h ^0 



\T{(vh)(jJ>h){Vh)\ 

\f>h\l\Vh\l 



But 



T[(v h )(ii h )(r) h ) = I D^h ■ Dr] h dx / [(cof D 2 v h )Dfi h ] ■ Dr] h dx 

Jn v Jn 

= [[(I- -(coiD 2 v h ))Dfx h } ■ Drjh dx, 
Jn v 

where / denote the n x n identity matrix. Hence 

'^y^ 1 =/[(/- i(cof &v h ))DJ*r) ■ D^- dx. 

Next, we note that for fixed Vh G X h , we can define a bilinear form on V h by the 
formula 



(P>?) 

Then since 



/ [(/ - -(cof D 2 v h ))Dp] ■ Dq dx. 
Jn v 



(p.?) = ^((p + q,p + q) - (p-q,p-q)), 



21 



we obtain 
\\T{(v h )\U 



< 



sup 



Vh 



\Vh\l \Vh\l 



[(I 



+ 



-cof D 2 v h )D( ^ 



- + - 


Vh n 


■D{^\ + 


1 


Vh\l 


\Hh\l 


Vh 


}-D( 


(J>h Vh 


\Vh\i 


\Vh\i \Vh\ 



dx 



Hh 


Vh 


) 


\V>h\\ 


\Vh\l 





a. 



Using ( |2.15[ ), we have 
(1 



v ' 



j j nrri 

w\l < [(I- -(cof D 2 v h ))Dw] ■ Dw dx < (1 )\w\\,w G #o(ft). 

Jo, u v 



Since i/ = (M+m)/2, 1-M/v = —(M—m)/(M+m) and 1-m/V = (M-m)/(M+m), 
we conclude that a < (M — m)/ (M + m). □ 

We can now prove the following lemma 

Lemma 2.12. The mapping T\ is a strict contraction in X h with contraction constant 
(M - m)/(M + m) for v = (M + m)/2. 

Proof. Let Vh and Wh G X h . Then, using the mean value theorem 



|Ti(iy h ) -Ti(v fc )|| 



< 



T[(v h + t(w h - v h ))(w h - v h ) dt\\ 
\T[{v h + t(w h - v h ))(w h - v h )\\dt. 



Since w h — v h e Vq 1 and + — t^) G X h , t G [0, 1], we obtain by Lemma 



2.11 



|Ti(w h ) - Ti(v h )\\ < 



1 M -m 
M + m 



w h - v h \i dt 



M - m, 
M + m 



□ 



Remark 2.13. For the operator T% to be a strict contraction, it is enough to have 
v sufficiently large. From the proofs of Lemmas 2.11 and \2.12\ we need to have max 
(\l—m/v\, \l — M/v\) < 1. This is equivalent to v > max (m, M/2) or M/2 < v < m. 

Remark 2.14. By the inverse inequality, we have 



C 5 h 2 \w h \l < \\w h \\l < C 6 \w h \l,w h G Vq. 



It follows that 

(C 5 h 2 --)\w h 
v 



r r | TTt 

\< w 2 h dx+ [-(cof D 2 v h )Dw h )-Dw h dx < (C 6 )\w h \\,w h G V h . 

Jo. Jn v v 



As in the proofs of Lemmas 



2.11 



and 



2.12 



we conclude that for max QC^h 



2 M I 



|C 6 - 



™|) < 1, the mapping T 2 : X h — )• (Vq)' defined by 



(2.35) < T 2 (v h ), fa) = [ v h i> h dx + - I (det D 2 v h - /) ^ dx, v h G X h , ^ G V h . 

Jn v Jn 
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is a strict contraction. But \C^h 2 — — | < 1 is equivalent to (i> < M / {C§h 2 )) or 
(v > M/(Cz,h 2 ) and h < 1/yfCl). And similarly, \Ce — — \ < 1 is equivalent to 
(y < m/(Co)) or (y > m/{C§) and 1 < 1/Cq). We conclude that T 2 is a strict 
contraction for v sufficiently small or for h sufficiently small and v in an appropriate 
range. 

We can now claim our main result, which is the convergence to Uh of the sequence 
defined by u k+ i G V h , u k+ i = g h on dfl and 



(2.36) v I Duk+i ■ Diph dx 



Du k ■ Dtp h dx+ (det D 2 u k - f) ijj h dx, ip h G V i 



h 

o • 



Theorem 2.15. Let Q be convex with a Lipschitz continuous boundary and assume 
that the spaces V h = S r d (l~) have the optimal approximation property (2.5) and satisfy 
the inverse estimates ( |2.7 ). The sequence defined by (2.36) converges to the unique 
strictly convex solution of (2.10) for any initial guess Uq in X h and a suitable 
v > with a linear convergence rate. 

Proof. The proof parallels Theorem 5.4 in [19]. Let us assume first that u k eX h . We 
have using ( |2.10 ), or equivalently det-D 2 -^ = / in (Vq)', 

I D(u k+1 - u h ) ■ Dip h dx = I D(u k - u h ) ■ Dip h dx + - f det D 2 u k ip h dx 
Jn Jn v Jn 

det D 2 Uh iph dx 



(ri(u fe ) - Ti(u h ),ij) h ). 



M -m, 



\u k - u h \i\u k+1 - u h \i, 



Taking iph = u k+ i — Uh, we obtain 

\uk+i -Uh\l < \\T\{u k ) -Ti(u h )\\ \u k +i —Uh\i < , r 

M + m 

where for simplicity, we assume that the finite dimensional V$ is equipped with the 
| \i norm of Hq(Q). We conclude that 



\u k+ i - u h \i < 



M -m, 



\Uk - u h U. 



M + m 

This also shows that if u k G X h , then u k+ \ G X h and concludes the proof. 



□ 



Remark 2.16. It follows from the above result and Remark \2.1J\ that for a suitable 
initial guess and a suitable u > 0, the sequence defined by 

(2.37) v / u k+1 tp h dx = v l u k ip h dx+ / (det D 2 u k — /) ip h dx,iph G Vq, 
Jn Jn Jn 

with u k+ \ G V h ,u k+ i = gh ondVt also converges to the unique strictly convex solution 
of 2.10. Obviously the convergence properties of (2.36) and (2.37) depend on the 
contraction constants ofT\ andT2 respectively. Thus (2.36) is more robust than (2.37) 
since the choice of v for (2.36) is less dependent on the discretization parameter h. 
As suggested in [26] in the context of monotone schemes, the use of the Laplacian 
preconditioner results in a more efficient algorithm. 



3. Numerical results 
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We give a brief review of the spline element method and give numerical results for both 
the latter and standard finite difference methods. The choice of the latter emphasizes 
the broad range of applications of the methods we proposed. The numerical results 
with finite difference methods complete the results with the spline element method 
in the degenerate case / > in Q and in three dimension. In that case, essentially at 
the discrete level / > Cq > for a constant Cq as only non zero values are used. In 
general we did not try to choose the value of v that would give the smallest number 
of iterations except in Tables [4] and [5] where we compare the performance of the two 
methods. 



3.1. Spline element discretization. We refer to [U El El El EEl E] for a description 
of the spline element method. We describe the method for linear problems and recall 



that the problems (1.3) are linear problems. Let u G V = H™(p,),m > 1 solve a 



variational problem a(u,v) = f(v) with the conditions of the Lax-Milgram lemma 



satisfied. Take Vh as the spline space S r d (T) of smoothness r and degree d, (2.4). For 
r = and d = 1 we have the space of piecewise linear continuous functions. 

First, start with a representation of a piecewise discontinuous polynomial as a vector 
in M. N , for some integer N > 0. Then express boundary conditions and constraints 
including global continuity or smoothness conditions as linear relations. In our work, 
we use the Bernstein basis representation, [U E] which is very convenient to express 
smoothness conditions and very popular in computer aided geometric design. Hence 
the term "spline" in the name of the method. We can therefore identify the space Vh 
with {c G R , Rc = G} for some integer N, matrix R and vector G. The discrete 
problem consists in finding c G Vh, c T Kd = F T d for all d G Vh for a suitable stiffness 
matrix K and a load vector F. Introducing a Lagrange multiplier A, the functional 

K(c)d - L T d + X T Rd, 

vanishes identically on Vh- The stronger condition 

K(c) + X T R = L T , 

along with the side condition Rc = G are the discrete equations to be solved. We are 
lead to saddle point problems 



K R T \(c 

r o Ma 



F 
G 



The ellipticity assures uniqueness of the component c and the saddle point problems 
are solved by a version of the augmented Lagrangian algorithm 

(3.1) (K + -R T R)c^ = K T c^ + —R T G, 1 = 1,2,... 

The convergence properties of the iterative method were given in [2]. Extensive 
implementation details can be found in [U E] ■ 
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h 




L 2 norm 


rate 


H 1 norm 


rate 


H 2 norm 


rate 


1/2 1 


236 


4.1569 l(T e 




6.5142 10" 5 




1.9364 10" 3 




1/2 2 


233 


1.1504 10" 7 


5.17 


2.3915 10" e 


4.77 


1.3444 10" 4 


3.85 


1/2 3 


233 


3.2406 10" 9 


5.15 


8.4120 10" s 


4.83 


8.9366 10" b 


3.92 


1/2 4 


233 


4.5857 10" 1U 


2.82 


4.7246 10~ y 


4.15 


6.0706 10" y 


3.88 



Table 1. Time marching method for Test 1, Sg, v = 50 



d 


n it 


L 2 norm 


H 1 norm 


H 2 norm 


3 


1 


1.2338 10~ 2 


7.6984 10~ 2 


4.4411 10" 1 


4 


270 


1.6289 10" 3 


1.4719 10~ 2 


1.3983 10" 1 


5 


135 


1.5333 lO- 3 


8.7312 lO" 3 


6.0412 10" 2 


6 


424 


1.2491 10" 4 


9.7458 10" 4 


1.0473 10" 2 


Rate 




0.18 O^" 1 


4.57 0.25 d 


60.85 0.3 d+1 



Table 2. Time marching method for Test 2 (3D) on Zi, v = 50 



d 


n it 


L 2 norm 


H 1 norm 


H 2 norm 


3 


1 


3.1739 10" 3 


2.3005 10" 2 


2.4496 lO" 1 


4 


651 


3.2385 10" 4 


3.5599 10" 3 


5.2262 10" 2 


5 


744 


2.2730 10" 5 


3.8977 10" 4 


8.8978 10" 3 


6 


652 


1.1956 10" b 


2.2056 10" 5 


6.0437 10" 4 


Rate 




0.72 0.072 d " 1 


29.44 0.1 d 


861.43 0.14 d+1 



Table 3. Time marching method for Test 2 (3D) on 7~ 2 , v = 50 

3.2. Numerical results with the spline element method. Unless otherwise indi- 
cated, all numerical simulations below are for r = 1 and the domain is [0, l] n , n = 2,3. 
For n = 2, the computational domain is the unit square [0, l] 2 which is first divided 
into squares of side length h. Then each square is divided into two triangles by the 
diagonal with negative slope. For n = 3, the initial tetrahedral partition 71 consists 
in six tetrahedra. Each tetrahedron is then uniformly refined into 8 subtetrahedra 
forming 7^. In the tables, n it denotes the number of iterations. We refer to [TJ [6] for 
implementation details of the method. 

We use test functions suggested in [THJ El ED] • 

Test 1: u(x,y) = e^ 2+2/2 ^ 2 so that f(x,y) = (1 + x 2 + y 2 )e^ x2+y2 ^ and g(x,y) = 

Test 2: u(x,y,z) = e ^ 2 +^ 2 )/ 3 so that f(x,y, z) = 8/8l(3 + 2(x 2 + y 2 + z 2 ^ 2 ^ 2 ^ 
and g(x,y,z) = e^ x +y +z ^ 3 on dVL. 

Barring roundoff errors, the methods introduced in this paper capture smooth solu- 
tions. For the two dimensional test function, Test 1, we give numerical results for 
successive refinements and for the three dimensional test function, we give numerical 
results for increasing values of the degree d on two successive refinements. 
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Figure 2. Time marching, Test 3, concave solution: h = l/2 A ,d = 
5,z/ = 50. 

In the context of approximations by finite dimensional spaces, many methods pro- 
posed, [HJ IH 120] , fail to fully capture the convexity of the solution on the test case 

Test 3: g(x,y) = and f(x,y) = 1. 

The pseudo transient method with L = A for plane problems enforces element by 
element Au > and for a non degenerate problem det D 2 u = f > 0. This implies 
that the numerical solution has Hessian positive definite element by element, which 
when combined with C l continuity gives numerical convexity [30], Lemma 1. The 
numerical results are given in Figure [4] with a plot of the graph of the solution along 
the line y = x. 

For the same test case, there is a concave solution. The concavity property of the 
concave solution obtained with the time marching method are better than the one 
obtained by the vanishing moment methodology, [20]. This is illustrated in Figure [2] 



We now discuss how the two methods compare. First, we are solving the same discrete 

Second, we noticed that the smaller 
Thus for a smooth solution, the correct 
value of v to take in the pseudo transient method is v = which is exactly Newton's 



equations (2.10) by different iterative methods 
v, the smaller the number of iterations 
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h 


V 


n it 


time 


L 2 norm 


rate 


1/2 1 





6 


3.032810 +u 


2.195410" 2 




1/2 2 





5 


8.136510 +u 


3.609710" a 


2.60 


1/2 3 





6 


3.823010+ 1 


1.068510 -3 


1.76 


1/2 4 


3 


56 


1.597910+ 3 


3.766610" 4 


1.50 



Table 4. Pseudo-transient method with L = A, Test 4 r = 1, d — 3 



ft, 


V 




time 


L 2 norm 


rate 


1/2 1 


2 


35 


6.2191 10+° 


2.072110- 2 




1/2 2 


2 


89 


6.0553 10 +1 


1.857910" 3 


3.48 


1/2 3 


4.5 


64 


1.6849 10+ 2 


5.043810" 4 


1.88 


1/2 4 


11.5 


151 


1.703810+ 3 


2.113210" 4 


1.25 



Table 5. Time marching method with L — A, Test 4 r = 1, d — 3 



method. In fact, Newton's method has been shown to have a quadratic convergence 
rate while the pseudo transient methods and time marching methods are shown in 
Theorems 2.9| and 2.15| to have a linear convergence rate. Moreover the numerical 



errors of Tables [TJ 2] and [3] are essentially the ones obtained with Newton's method 
as expected. We compare the performance of the methods on a non-smooth solution 
with known solution. 

Test 4: u(x,y) = —a/2 — x 2 — y 2 with corresponding / and g. 

The time listed is in seconds and obtained on an imac running Mac OS 10.6.8 with a 
2.4 Ghz intel core 2 duo and 4 GB of SDRAM memory. While for small values of h the 
time marching method appears to take significantly more time, it is also significantly 
more accurate. For h = 1/2 4 the time took by the two methods is almost the same 
with the time marching method giving a more accurate solution. 

Note that the convergence analysis in Holder spaces require the domain to be smooth 
while the convergence analysis in Sobolev spaces only assumes the existence of a 
smooth solution with the domain allowed to have Lipschitz continuous boundary. We 
conclude this section with a test problem on a non square domain to contrast with 
results that can be obtained with the finite difference methods. 

Test 5: we consider the unit circle discretized with a Delanauy triangulation with 824 
triangles and u(x,y) = x 2 + y 2 — 1 which vanishes on the boundary, Figure K3j 



3.3. Numerical results with finite difference methods. We used a standard 



finite difference discretization of (1.3). The results are given on Tables [6j [7] and 
Figure |4j Numerical errors are in the maximum norm. The method also behaves 
well, in the non degenerate case, for the three dimensional analogues of these cases 
(also used in [23J). 



For u{x,y } z) = — y 3 — x 2 — y 2 — z 2 , it was necessary to stabilize the method using 
( 1.4[ ) with m large. With m = 25, v = 50, we obtained the same results with the 
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Figure 3. u(x,y) = x 2 + y 2 — 1 on a non square domain with pseudo 
transient u = 0,r = l,d = 3 

h 



V 


1/2 2 


1/2 3 


1/2 4 


1/2 5 1/2 6 


50 


3.9093 10- 3 


1.0340 10~ 3 


2.6643 10~ 4 


6.6964 10~ 5 1.6781 10~ 5 



Table 6. Smooth solution u(x,y) = e^ x2+y ' 2 ^ 2 
h 



V 


1/2 4 


1/2 5 


1/2 6 


1/2 7 


1/2 8 


15 


2.211310-* 


1.5945 10- 1 


3.8944 10 u 


6.8634 10 +i 


1.1695 10+ 3 


150 


2.2113 10~ 2 


1.6920 10~ 2 


1.2440 10" 2 


3.3583 10~ 2 


1.1316 10° 


250 


2.2113 10~ 2 


1.6920 10~ 2 


1.2440 10~ 2 


8.9702 10- 3 


2.1435 10" 1 


600 










1.6745 10~ 2 


800 










6.4054 10" 3 



Table 7. Non smooth solution (not in H 2 (Q)) u(x,y) = — a/2 — x 2 — y 2 



three dimensional iterative method introduced in [I], 

Au k+1 = ({Au k ) s + 9{f -detD 2 u k )^. 

The errors in the maximum norm were given by 3.0976 10~ 3 , 1.0432 10 -3 , 1.4169 10 -3 , 
1.3766 10- 3 , 1.1017 10- 3 , 8.3671 10" 4 , 3.6635 10~ 5 for h = l/2 k ,k = 2,...,8 respec- 
tively. 

The convergence rate for the smooth solution is 2, and 0.50 for the two dimensional 
non smooth solution, and 0.25 for the three dimensional non smooth solution. For 
the last one, we took only in account errors for h = l/2 fc , k = 4,5, 8. 
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Figure 4. No known exact solution (left), f(x, y) = 1, g(x, y) = 0, h = 
l/2 5 ,z/ = 15 and u(x,y) — \x — 1/2 j with g(x,y) — \x — 1/2 1 and 
f(x,y) = 0, h = l/2 2 ,u = 5 
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